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ABSTRACT 

We present a numerical model in which a cold pair plasma is ejected with relativistic 
speed through a polar cap region and flows almost radially outside the light cylinder. 
Stationary axisymmetric structures of electromagnetic fields and plasma flows are self- 
consistently calculated. In our model, motions of positively and negatively charged 
particles are assumed to be determined by electromagnetic forces and inertial terms, 
without pair creation and annihilation or radiation loss. The global electromagnetic 
fields are calculated by the Maxwell's equations for the plasma density and velocity, 
without using ideal MHD condition. Numerical result demonstrates the acceleration 
and deceleration of plasma due to parallel component of the electric fields. Numerical 
model is successfully constructed for weak magnetic fields or highly relativistic fluid 
velocity, i.e, kinetic energy dominated outflow. It is found that appropriate choices of 
boundary conditions and plasma injection model at the polar cap should be explored 
in order to extend present method to more realistic pulsar magnetosphere, in which 
the Poynting flux is dominated. 
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1 INTRODUCTION 

A global structure of pulsar magnetosphere is one of key issues to understanding energy outflow to the exterior. The numerical 
model has been successfully developed in the past decade, although the basic equatio n was already derived i n the early da ys 
of pulsar theory. Extens ive reviews are available in some books (e.g, iMichell (|l99lh : iBeskin et all (|l993l ): iMestell ^^). 
Contopoulos et al. 1 19991 ) calculated the stationary axially symmetric magnetosphere based on the force- free approximation. 
They for the first time showed a solution with dipole magnetic field lines near a neutron star, which smoothly pass through 
the light cylinder to the wind region at infinity. In the model, there is a current sheet fiowing on the separatrix and equator 
outside the light cylinder. The magne tosphere model is subs equently explored in detail by several authors; some physical 
properties re presented by the solution JOgura fc Koiima 20031 ). the Y-point singularity between open and clo sed field lines on 
the e quator ([Uzdenskvlbood : iTimokhinl bood). the electromagnetic luminosity in high numerical resolution l|Gruzinovll2005l . 
20061 ). Numerical construction of the magnetos phere around an a l igned rotator is also performed using time-dependent codes 
with the force- free and MHD approximations IIKomissarov 20061 : McKinnev 20061 ) . A stationary state, which is very similar 
to the solution given by Contopoulos et al. ( 19991 ), is ob tained with certa i n initial and boundary conditions. Th e approach is 
extended to an oblique rotator by 3D simulation codes ( Spitkovskv 20061 : Kalapotharakos &: ContopoulosI 2008h . 

An ideal MHD condition E — B x v/c is used to determine the electric field in these calculations irrespective of the 
numerical methods. Consequently, two electromagnetic field vectors are always orthogonal B ■ E = 0, and the parallel com- 
ponent of E along the p lasma motion vanishes e verywhere. This condition holds if the plasma density much exceeds the 
Goldreich- Julian density ( Goldreich fc Julian 19691 ). The global structure based on the force- free and MHD approximations is 
obviously a good step to understanding the whole magnetosphere. However, it is important to study how and where the con- 
dition breaks down, and how this changes the electromagnetic field structure and plasma behavior. An alternative approach, 
in which the ideal MHD condition is relaxed, is necessary to address these problems. Breakdown of ideal MHD condition in 
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the pulsar magnetosphere is qualitatively pointed out in the literature, e.g, (jMestel fc Shibatalll994l : iGoodwin et al.ll2004l l. 
It is our purpose to further study the problem by actual modeling. It is necessary to determine the electric fields from the 
distribution of the charge density, since E = B x v/cis no longer used. The electric acceleration or deceleration of fluids will 
be allowed elsewhere, since v ■ E ^ 0. The location may be important to the observation. 

In this paper, we present an approach based on a two-fluid plasma consisting of positively and negatively charged particles. 
In this approach, the electromagnetic fields are modeled by Maxwell's equations with a plasma source. Resultant electric fields 
are not in general perpendicular to magnetic fields B ■ E ^ 0, i.e, breakdown of ideal MHD condition. The stream lines of the 
plasma flows are not a priori assumed to coincide with the magnetic field lines. The plasma flows are determined by equation 
of motion under the electromagnetic forces. In this paper, we consider a simple plasma model, a cold dissipationless plasma, in 
which thermal pressure, pair creation or annihilation and radiation loss are neglected. If we furth er neglect the inertial terms 
in equations of motion, then we have the force- free and MHD conditions ([Goodwin et al.ll2004l ). We here keep the inertial 
terms, in order to model the magnetosphere taken into account one of non-ideal MHD effects. We also assume stationary 
axisymmetric states in the electromagnetic fields and the plasma fiows. However, the azimuthal component of a vector is not 
zero in general. For example, a toroidal magnetic field B^ can be produced by poloidal currents jp — {jr,je), in the spherical 
coordinates (r, 6, This approach can be naturally extended to include other physical mechanisms in the future. This paper 
is, therefore, the first step towards including more physical processes. 

This paper is organized as follows. In section 2, we discuss our numerical method to solve electromagnetic fields and fluid 
streams in two-dimensional meridian plane. The relevant boundary conditions are also given. A model of plasma injection is 
given at the inner boundary, a polar cap region. The angular dependence of the injection rates is different between positively 
and negatively charged particles, but the total is the same. The total current from the polar cap is therefore zero. This provides 
a model to construct magnetosphere with charge-separated plasma flow. This choice of the injection model is not unique. In 
section 3, we show our results of the global structure. Section 4 contains our conclusions. 



2 ASSUMPTIONS AND EQUATIONS 

2.1 Electromagnetic fields and plasma fiows 

Axially symmetric electromagnetic flelds in the stationary state are expressed by three functions ${r,6),G{r,9),S{r,8) as 
E = -V-1', (1) 
B = IvG X e-^ + |e-<^, (2) 

where we use spherical coordinate (r, 0, </)) and R = rsind. It is convenient to use the non-corotational potential 'I' = $ — JlG/c, 
where fl is angular velocity of a central star. Maxwell equations with charge density pe, and poloidal and toroidal components 
of current {jp,j,p) are given by 

VG^-^RU, (3) 

1 ^ ^ - 47r ^ , , ^ 

_VS X = —Jp, (4) 

^2r , ( ^R .\ 2n f dG ^dG\ 

where operators V and in spherical coordinate are given by 



^ sine d 



V 



2 I d f 2 d \ . 1 d f . „d 



r.2 



We adopt a treatment in which the plasma is modeled as a two-component fluid. Each component, consisting of positively 
or negatively charged particles, is described by a number density n± and velocity v± — v±p + v^e^. Note that the proper 
dens ity n+ is related with the lab-frame density n± by n% = n±/7±, where 7± is a Lorentz factor 7^ = (1 — {v±/cf')~^^^ 



(e.g. iGoodwin et al.l (120041 )). We assume that the positive particle has mass m and charge q, while the negative one has mass 



m and charge —q. The charge density and electric current are given in terms of n± and v± as 

= (7(n+-n_), (8) 
J = q{n+v+ -n-v-). (9) 



Continuity equation for each component in the stationary axisymmetric conditions is 
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= V • {n±v±) = V ■ {n±v±p). 

The poloidal velocity components v±p are satisfied by introducing a stream function F± (r, 6) as 



n±v±p 



R 



VF± X 



From the definition, the number density is given by 
^ |VF±| 

From eqs.® and pijl . the current function 5* in eq.Q can be solved as 



(10) 



(11) 



(12) 



(13) 



The electromagnetic force is dominant so that collision, thermal pressure and gravity are ignored. The interaction between 
two-component fluids is assumed only through the global electromagnetic fields. The equation of motion for each component 
with mass m and charge ±g in the stationary state is given by 



V) 7± 



v± 



E 



v± 



X B 



(14) 



By adding and subtra cting equations (1141) for two components, we hav e an equation of one-fluid bulk motion and a generalized 
Ohm's law. See e.g, ( Melatos fc Melrose 19961 : Goodwin et al. 2004 ) for the detailed discussion. We do not follow such a 
treatment, but rather solve ea. (|14|l for each component. Using the identity (v ■ V)jv = (V x 7«) x « -I- V7C^, we find 
two conserved quantities along each stream line, corresponding to axially symmetric and stationary conditions. They are 
generalized angular momentum J ± and Bernou lli integral K± , which are obtained by the azimuthal component of eq. ([14} and 
a scalar product of v and eq. mi) (|Mestellll999h . Their explicit expressions are given by 

q 



J± 



'y±v±d,R ± — G, 
mc 



(15) 
(16) 



These quantities depend on the stream functions F± only, and the spatial distributions are therefore determined by F± which 
is specified at the injection point in our model. The convenient form for the third component of ea. (|14p is the azimuthal 
component of a cross product of v and eQ. (|14|l . This means a perpendicular component to the stream lines, which is given by 



2?F± = Vln( ^ 

,7± 



7± V c^R / mc 7± 



where J'± and ifj- are derivatives of J± and K± with respect to F±. Using VJ± 
ea. (|17|l can be written in an alternative form 



VF± = 



Vln 



n± 

7± 



7±(i^±r + f±e) 



mc 7± 



(17) 



J±VF±, VK± = K'^VF± and eq.Hl, 



(18) 



The stream function can not be defined in corotating region, where poloidal components of the velocity vanish, and the 
expression (|ll|l is no longer used. Instead, the charge density pe* and current j, are given in terms of the corotating condition 
$ — QG/c, j, — p^fQ.Re4,. From eq.([5l, the corotating charge density is given by 

Three velocity components, Lorentz factor and number density are determined by ea. (|12|l and two integrals (|15|) . (|16p . 
if four functions G, ^' and F± are known. The charge density (|8} and current ((9)1 are calculated from these fluid quantities 
of both species. In the corotating region, they are given by corotating charge density and current. Irrespective of the spatial 
region, the source terms of partial differential equations for G, ^ and F± depend on themselves in a non-linear manner. Some 
iterative methods are needed to self-consistently solve a set of eqs.([3l,(l5} and (|18p . There is no established method so far to 
solve nonlinearly coupled equations, so that our numerical procedure is rather primitive. Initial guess for these functions is 
assumed, say G'"', and f!^\ Using these functions, the source terms are calculated, and a new set of functions G'-^^ 



and F^"' are solved from these source terms with appropriate boundary conditions. The procedure is repeated until the 
convergence, say, |G^"+^' - G'"'|, |*("+^^ - ^f'"'!, \F^+^^ - < e, where e is a smaU number. The iteration scheme may 

not necessarily lead to a convergent solution, since there is no mathematical proof. 

In order to e xamine our n umerical scheme, we have performed a test for the split-monopole case, for which an analytic 
solution is known (|Michelll 19731 ). The non-corotational electric potential in the solution is zero everywhere, so that the condition 
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Table 1. Convergence test 



Model 


grid 


II^GIIo 


P^+llo 




ll<5G||iv 




¥p+\\n 




-\\n 


Ai 


150 X 50 


1.0 


0.15 


0.15 


1.4 X 10" 


-3 


1.3 X 10-2 


1.1 X 


10-2 


Aa 


300 X 100 


1.0 


0.15 


0.15 


4.1 X 10- 


-4 


7.4 X 10-3 


5.6 X 


10-3 


Bi 


150 X 50 


1.0 


0.25 


0.25 


2.5 X 10" 


-3 


1.6 X 10-2 


1.2 X 


10-2 


Ba 


300 X 100 


1.0 


0.25 


0.25 


9.1 X 10" 


-4 


7.8 X 10-3 


5.2 X 


10-3 



^ = is used and a reduced system of G and F± is checked. These functions are numerically solved by a finite difference method 
with appropriate boundary conditions in the upper half plane. Results for the convergence to the solution are given in Table 
[T] Two types of initial trial functions and two different grid numbers are used. Deviation from the analytic solution is shown 
by a norm ||<5/|1„, which is evaluated at all grid points as = ^j) -,/* (n, ' / [E(/* . 

where /* is the analytic solution and is numerical result after n iterations. We have repeated until the relative error 
e = 1 X IQ-'^ in this test problem. We have started from G'^' = 0, so that ||<5G||o = 1. The initial choice of G^"-* is not so 
important, since the numerical solution approaches the analytic one at the first step. On the other hand, the choice of initial 
guess for F± is important, ft is not easy to set large deviation at the initial step, since the function F± should be monotonic. 
If there is a maximum or minimum, where |VJ±| = inside the numerical domain, the flow vanishes nvp = 0. This causes 
a numerical difficulty at that point. From the monotonic nature consistent with the boundary conditions, the initial norm 
||5J±||o can not be large. Table[T]shows that the numerical solutions successfully converge on the analytic ones within certain 
errors. Convergence factor e does not exactly correspond to deviation from true solution, but gives an estimate. The true 
solution is not known in most problems, and the deviation can not be calculated. The convergence factor e can be regarded 
as error estimate. 

2.2 Boundary conditions 

We assume that the axially symmetry around a polar axis and the reflection symmetry across an equator. The numerical 
domain in the spherical coordinate {r,6) is ro < r < ri, < S < 7r/2. The inner and outer radii in our calculations are 
ro = tl /5 and r\ — Sr^ , where = c/il is the distance to the light cylinder. Figure 1 schematically represents the numerical 
domain. The functions at the inner boundary ro are closely related with plasma injection model, which is separately discussed 
in the next subsection. We here discuss the boundary conditions at the axis, equator and outer radius. 

We solve the magnetic flux function G in the upper half plane between ro and ri , which is the region enclosed by a curve 
PBRQP in Fig.l. Poloidal magnetic field at the inner boundary tq, i.e, on PB is dipole, so that we impose the condition, 
G = ^s\v? 9/ro, where /i is the magnetic dipole moment. The plasma is injected through a polar cap region Q < 9 < 9o aX 
ro, PA in Fig.l. A curve AL represents the last open magnetic field line. All the field lines originated from a point with 
9 < 9o &t ro extend to infinity, whereas those from 9 > 9o are closed. The point L is (rL, ■7r/2). The last open line is given 
by Go ~ fJ. s'm^ 9o/ro. For purely dipolar field, the polar cap region and critical field line are given by sin2 6o — ro/vL and 
Go = ft/rL- These values in our numerical model are not known a priori, but are determined simultaneously with the global 
structures. The boundary condition on the equator 9 = n/2 is dG/d9 = inside the light cylinder, i.e, on BL in Fig.l, while 
outside it G = Go. This condition on the equator means Br = Q inside the light cylinder, but Bg =Q outside it. The boundary 
condition at the outer radius on QR is continuous, dG/dr = at ri. This condition means that the poloidal magnetic field 
becomes radial, since Be = 0. On the polar axis PQ, we impose the regularity condition which is given by G oc (rsinS)2 for 
9 ^ Q, i.e. Be = on the axis. 

Next, we consider the boundary conditions for the stream function F±, which is defined outside of the corotation region, 
i.e, a region enclosed by a curve PALRQP in Fig.l. We here assume that irrespective of the fiuid species, the last stream 
line coincides with the last open line of the poloidal magnetic field. The function F± should continuously approach a constant 
F±c on the boundary with the corotation region AL. Outside the light cylinder, the boundary condition on the equator LR 
is F± = F±c for rL < r < ri at = 7r/2, since there is no flow across the equator due to the refection symmetry. Outer 
boundary condition at ri is dF±/dr — 0. This condition also means that the flow becomes radial since nvg — {dF/dr)/R = 0. 
The regularity condition on the polar axis PQ is F± oc (r sin 0)2 for 9^0. 

Finally, we consider the boundary conditions for 'i', which is non-corotating part of the electric potential. We solve it 
only outside the corotating region, the region enclosed by a curve PALRQP in Fig.l, since = in the corotating region. As 
the boundary condition of ^P, the function contirmously becomes zero, * = 0, toward the boundary AL with the corotation 
region. Outside the light cylinder, the condition on the equator LR is assumed as = at S = 7r/2 for tl < r < ri. Outer 
boundary condition at ri is also assumed as = 0. One might think the boundary condition on the polar axis PQ is 'i' — 0, 
for which ideal MHD condition B ■ E — is satisfied due to E,. — Bg — 0. We found that the Dirichlet condition is too severe 
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O B L R 

Figure 1. Meridian region of numerical calculation. The region enclosed by a curve ABLA is a corotating region. Plasma is injected 
through the polar region PA and goes out through the outer radius QR. A line OQ is polar axis, and OR is an equator. 



to lead to any numerical solution. We use the regularity condition of 5', that is, /dO = 0. The Neumann condition is less 
severe, and allows the numerical solution. This means Eg = Q on the axis. The ideal MHD condition is likely to be broken 
near the axis, B ■ E ^ 0. This condition is quite diflterent from usual MHD treatment. 



2.3 Injection model 



In our model, plasma is assumed to flow through a small polar region, < 9 < 9o, ai the inner boundary ro, represented 
by PA in Fig.l. If the stream lines completely agree with the magnetic dipolar field lines, then the stream function is given 
by F± oc sin^ 6/r near the surface ro. We assume that the function for each particle type slightly deviates from the dipolar 
configuration near the polar cap region. In our injection model, the stream function in the range of < < is given by 

A sin^ ( 



F±iro,e) = 



ro 



1=F« 1 



/ sin 61 y 
V sin 6^0 / 



and the current function S is calculated as 

/ N I M r7 ' ^ 

1 



SttqAq . 2 
sm ( 



cro 



sm U 
sin 8o 



(20) 



(21) 



The constant a determines the deviation from the dipolar field. The poloidal current completely vanishes in the limit of 
a — 0, where the positively charged particles and negatively charged particles move along the common stream lines. The scale 
factor A is chosen as AnqXa = /ifi, so that the current function can be written as S ~ —2^lGd/c near the polar region, where 
Gd = fisin^O/r. The current function corresponds to the split-monopole solution near the polar region (Michel 1973). Our 
current function S smoothly goes to zero at the edge of polar cap, 9 = 9o. T his property is different from that of the force-free 



and ideal MHD approximations in which the function has a discontinuity ( Contopoulos et al. 19991 : Ogura fc Koii: 
Gruzinovll2005l ) . That is, there is a current sheet. 



20031 : 



The current density as a function of angle is calculated as jr — —4:fiflcos9[l — 2(sinS/ sin6'o)^]/ro at the polar cap. The 
electric current is negative for < sin 61 < sin&o/\/2 and positive for sin0o/\/2 < sm9 < sin 6*0. The total current ejected 
through the polar cap region is zero, since S{ro,9o) = 0. The positive or negative current flow is produced from the charge- 
separated plasma, i.e, different number-density distribution between two components in our model. The injection flow at ro 
is calculated from ea. (|20|l as 



n±v±r 



2A cos t 
r'o 



l=Fa 1 



Y sin 6* y 
V sin 6*0 / 



We assume that the flow speed vq is relativistic, 70 — {I — {vo/c)'^)~^^^ 2> 1, and is independent of ( 
to be the same for each particle type. The flow direction through the small polar region cos 6 ~ 
velocity is Vr ~ c, so that the number density at ro is approximately given by 



2A 



n± = 



iTQ 1-2 



/ sin 6* y 
V sin 9q J 



(22) 

This property is assumed 
1 is almost radial and the 



(23) 



The charge density is given by 



AqXa 






erg 


V sm 9q J 





/ sin 9 \ 
V sin ^0 / 



(24) 
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where our choice of parameter AnqXa = fiQ is used. The typical value of ea. (l24|l /if2/(7rcro) = Bd^i/{2'Kc) is the Goldrich- Julian 
charge density for the field strength Bd of magnetic dipole. 

Near the polar cap region, the force-free condition is satisfied, so that the current function 5* and electric potential ^' 
depend on the magnetic flux function G. We adopt the following forms, Sp{Gd) and $p(Gd), as a function of dipolar flux 
function Gd as 

S,^-^^{l-bGd), (25) 
c 

= ^ + ^ [1 - ibGdf] . (26) 

where 6 — ro/{fism^ 9o). Equation (|25p is reduced to ea.(|21|) at ro, and the Poisson equation is approximately satisfied for $p 
(|26|) and the charge density (|24l) . The electric current in the force- free condition is generally given by 

j^,_B + p.cR-e,. (27) 

By the straightforward calculations, it is found that the poloidal components of ea. (l27p with the expressions Sp{Gd) and 
^p{Gd) are satisfied and toroidal component gives a small value « p^c x (ro/r^)^''^ ^ jV ~ PeC. We regard = and 
impose the ^-component of the fiuid velocity as 

v+,i, = — — Jto", v-,/, = — ^ -Koii, (28) 

where -Ro ~ ro sin 6. For this choice, total angular momentum of plasma fiow is m{n- + n+)-yoRo^ at the inner boundary. 

We here summarize the boundary conditions at ro. Equation (|20|) is used for F± with A-KqXa = fiQ,, and eq. (|26|) for $. 
Two integrals J± and K± are calculated at ro as a function of polar angle 6 from eqs. (|23p . (|26p . (|28|l . Gd and 70. The relations 
J±{F±) and K±{F±) are constructed by eliminating 6 in terms of ea. (|20|) . 



3 NUMERICAL RESULTS 



We use a finite difference method to solve a set of partial differential equations. The typical grid number in the spherical 
coordinate {r,9) is 300 x 100 for 0.2 < r/rL < 5 and < 6 < tt/2. The polar cap region at the inner boundary is covered 
by approximately 30 grid points. We have obtained the same result by changing the grid numbers as 150 x 50 or 450 x 150. 
The convergent factor is e ~ 10~^, and can not be improved so much by the grid refinement. We demonstrate numerically 
constructed magnetosphere. Parameters used in the numerical calculation are a deviation parameter a — 0.2, Lorentz factor 
70 = 10^ and q/i/ {mc?r\) = 10. The last dimensionless parameter is magnetic gyration frequency to angular velocity Q. = c/tl 
of a star. Numerical results of the plasma flows and electromagnetic fields depend on a single combination, 77 = qfi/ {mc^'yorj^), 
as far as 70 ^ 1, since the source terms for eqs.©,© and pSf) are scaled by it. Thus 77 is a key parameter to determining 
the global structure. See Appendix for the details. 

Figure [2] shows numerical solution of the magnetic function G. We also show that of dipole field Gd = p.siv?' 9/r for the 
comparison. The poloidal magnetic field is dipole near the inner radius ro/rL = 0.2. The field gradually deviates from the 
dipole, and becomes open outside the light cylinder. The field configuration is eventually radial near outer radius ri/r_L = 5. 
The numerical result provides the la st open field line as Go ~ 1.15^/rL. The c ritical value is Go = 1. 5927j/rL in magnetospher e 
fulled with rigidl y rotating plasm a (lMichellll973 . Il99ll : iMestel fc Prvcelll992h . and Go = 1.36/i/rL JContopoulos et al.lll999l '). 
Go = 1.27///ri, ( Gruzinov 2005h in a solution with the force-free approximation, and Go = 1.26fi/rL ( Komissarov 2006 ) in 
MHD simulation. Our value is smaller than that of other models, but is not so different. The total current is different from 
that of these models, so that there is no reason why the critical value should agree. 

Figure [3] shows the results of the stream functions for both species. The global structure of stream lines is almost the 
same as that of the magnetic field lines shown in Fig[2] although the numerical agreement is not so complete. Thus, fiows 
in meridian plane is almost parallel to the magnetic field lines. A difference between F+ and F- at large radius originates 
from the inner boundary condition at ro, where the fraction of negatively charged plasma is slightly larger in polar region 
« 0, but smaller for ~ Oq. The property extends to the outer radius. At large radius, the flow becomes is radial and the 
velocity is still relativistic, so that the number density decreases with the radius, n± ~ \dgF±\/{r'^ sin 9) oc r~^. However, the 
fraction (n+ — n_)/(n+ + n_) is still finite, and the charge separation remains. Our numerical model shows that negatively 
and positively charged regions are separated approximately by a curve with F± ~ 0.5(A sin^ ^o/j'o). 

We show numerical result of the electric potential. The contour of the non-corotating part = $ — ilG/c is shown in 
Fig. m There is a peak on the polar axis at r ~ 0.4rL. The function decreases toward the outer boundaries, where vl> = is 
imposed as the boundary condition at outer radius, and on the last open magnetic field. Numerical result shows the maximum 
value ^ = 0.85/ir2/(crL) at {r,0) — (0.4rL,0). This value is not small, since the maximum of corotating electric potential is 
Q.G/c ~ 1.15/-if2/(crL). Total electric potential $ = "I* -f QG/c is shown in Fig. [5] Overall structure is very different from the 
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Figure 2. Magnetic flux functions. Solid curves de- 
note the flux surfaces of numerical solution for Gri^/fi = 
0.1, 0.2, ■ • • , 1.4, and dotted curves those of dipole for 
GdfL/fJ- = 0.2, 0.4, ■ ■ ■ , 1.4 starting from the polar axis. 



Figure 3. Stream functions for negatively and positively 
charged fluids. Solid curves denote the flux surfaces of F_, 
and dotted curves those of in intervals of 0.1(Asin^ do/ro). 




magnetic flux function G or stream functions F± shown in Figs[2El The difference is clear at the polar region, whereas the 
agreement becomes better at high latitude region near the equator. Ideal MHD condition B ■ E = is not assumed in our 
model. The deviation is very large on the polar axis. This feature is closely related with the boundary conditions. As discussed 
in section 2, the boundary condition of 'i' is not Er — 0, but Eg = on the axis. This mathematical condition may allow 
the large value of 'I' on the axis. On the other hand, ^' almost remains zero near the equatorial region from the boundary 
condition. 

We discuss a consequence of non-ideal MHD field B ■ E ^ 0. Figure |6] shows contour of the Lorentz factor of positively 
charged particles normalized by initial one 70. At the injection boundary, 7+ is fixed as 7+ = 70 = 10^ for all polar cap 
angle, but there is a gradual increase toward the outer radius. The increase is remarkable at low 6, but 7+ is almost constant 
for the flow along the equator. The increase of 7+ is determined by the Bernoulli integral as A7+ — —qlS.^/{mc?), since 
= /S.K+ — A7+ -|- gA$/(mc^) along each flow line. Thus, large acceleration of positively charged particles toward the polar 
region can be understood, since available potential difference — A"l> is large as inferred from Fig. [5] The change of the Lorentz 
factor 7_ of negatively charged particles is opposite in the sign, and is given by A7_ — gA$/(mc^). They are therefore 
decelerated toward the polar region. 

Figure [7] demonstrates electromagnetic forces acting on a positively charged particle on the stream lines. Two vectors qE 
and q{v x _B)/c in the meridian plan are shown by arrows. The sum of these two forces causes a net acceleration. It is clear 
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that the vector E is not perpendicular to the flow hnes, and hence accelerates outwardly. The acceleration mechanism works 
better for the flow toward the polar region. This is another explanation for the increase of the Lorentz factor 7+ as shown in 
Fig. m This electric effect is opposite for negatively charged particles, which should be decelerated. 

In Fig. [S] we show the current function S. The poloidal current flows along a curve with a constant value of S. Figure [8] 
demonstrates a return current. That is, two distinct positions at ro are connected by a curve, say, the curve with 5 = -0.3p/ri. 
Such a global return current is generally produced due to the inertial term in our model. The integral J± in ea. (ll5|l is replaced 
by G in the limit of m = 0. The stream function F± is constant on a constant magnetic surface G, and the current function 
S = 47rg(_F+ — F-)/c is also constant. The global structure of S should be the same as that of G. Therefore, no loop of S is 
allowed in the limit of m = 0, since G is open field outside corotation region. In the context of generalized Ohm's law, the 
inertial term is a kind of resistivity. This term causes the dissipation of global current, and return current is produced in our 
model. 

We consider the effect of the current decay on the toroidal magnetic field, which is given by — S/R. Figure |9] shows the 
global structure. The function is zero at the polar axis, on the last open field line of G, and has a maximum at ^ ~ 60 /V^ of 
the polar cap region. The function decreases outwardly. Ratio to the poloidal component is important since the magnetic field 
strength both of poloidal and toroidal components decreases with radius. The ratio is small, \B^/Bp\ < (ro/r_t)'^''^/2 ~ 0.04 
at the inner boundary. We numerically estimated and found that \B^/Bp\ « 0.5 at (r_t,7r/4), and \B^/Bp\ « 1 at (4r_L,7r/4). 
Outside the light cylinder, the poloidal magnetic field is monopole-like as shown in Fig. [2] so that |_Bp| cx r~^. On the other 
hand, jS^j oc along the stream lines in the limit of m = 0. The slope of \B^\ slightly becomes steep due to the inertial 
term, but is not so steep as oc . In this way, toroidal component of the magnetic field is gradually important with radius, 
although the resistivity is involved in our model. 

The electromagnetic luminosity through a sphere at r is evaluated by radial component of the Poynting flux as 

Lom(r) = 2 / —(Ex B)r2TTr^ sinOde. (29) 
Jo 47r 

The luminosity of the plasma flow is a sum of both species as 

Lpis.sui3.{r) = 2 / mc^{"f+n+v+r + 'y-n-v-r)2nr'^ sm9d6. (30) 
Jo 

The energy conversion between two flows is possible through the Joule heating j ■ E, but the total Lem(T') + ipiasma(»') should 
be conserved. In Tabled numerical results are shown for different radii. We can check the conservation within a numerical 
error. The energy flux by plasma flow is always much larger than the electromagnetic one, and the conversion is very small 
in our model. The magnitude of the luminosities is almost fixed by the injection condition. We analytically evaluate these 
luminosities at ro using the inner boundary conditions, and find that Lcm = (/lilsin^ 9o/ro)^ /{3c) ~ 0.4/i^f2*/c^ and Lpiasma ~ 
87rmc^7oAsin^ 6o/ro ~ 2^^Q?' siv? 6o/{ar]ro) ~ llS/i^Sl^/c^, where the numerical values a — 0.2, rj = 0.1 and sin^ 60 = 1.15f2/c 
are used. In order to simulate the Poynting flux dominated case, it is necessary to increase the parameter 77. 



Magnetosphere with Two-fluid Plasma Flows 9 




Table 2. Luminosity through a sphere with radius r. 



1.5 0.18 115.57 

2.0 0.20 115.55 

2.5 0.21 115.54 

3.0 0.22 115.53 

3.5 0.23 115.53 

4.0 0.23 115.53 

4.5 0.23 115.53 



4 CONCLUSION 

We have numerically constructed a stationary axisymmetric model of magnetosphere with charge separated plasma outflow. 
The stream lines of pair plasma are determined by electromagnetic forces and inertial term. The massless limit corresponds to 
the force-free and ideal MHD approximations. The global structures of electromagnetic fields and plasma flows are calculated 
by taking into account the inertial term. In particular, the non-ideal MHD effects are studied. The electrical acceleration or 
deceleration region depending on the charge species appears. Poloidal current slightly dissipates. Numerical results depend on 
a single parameter rj = q^/ {mc?^ar\) as far as 70 ^ 1. The number of our model demonstrated in section 4 is 77 = 0.1, and is 
small when applying to the pulsar magnetosphere. Typical number 77 is estimated for electron-positron pair plasma, magnetic 
field Bs at the surface and spin period P as rj = 10'*(i3s/10^^G)(P/ls)~^(7o/10^)~^. Our present model is only applicable to 
highly relativistic injection (70 ^ 1) or weaker magnetic fields {Bs ^ 10^^ G). It will be necessary to scale up many orders 
of magnitude to 77 ~ 10* in order to apply the present method to more realistic cases. We in fact tried scaling up in our 
numerical calculations, but found that it is not straightforward. 

The difficulty and the limitation to smaller value of r) are closely related with boundary conditions and involved physics, 
as explained below. The Bernoulli integral (|16|) should satisfy a constraint K± =F q'^ /(rn(?) = 7 > 1. In our model, we specify 
K± by the injection condition which is fixed at the inner boundary. During the numerical iterations, the magnitude of the 
function ^ becomes very large in a certain region, where the condition K± =p q'^ /{mi?) > 1 is no longer satisfied. It is easily 
understood that this easily happens for large value of = g/i/ (mc^ 70 r|,), since the typical scale of is ji/r\. This gives 
a certain upper limit to the choice of 77. The actual estimate of the limit is somewhat complicated, since the potential ^ 
depends on the choice of boundary conditions, especially injection model at inner radius. It is therefore necessary to explore 
consistent boundary conditions or to include some physical process, in order to calculate the models with larger rj. Adjusting 
mechanism may be required at the boundary or within the numerical domain. For example, our present inner boundary is 
one-way, i.e, injection only at fixed rate. During numerical iterations, the charge density might numerically blow up elsewhere 
due to poor boundary condition. If the injection rate is able to be adjusted or plasma is absorbed through the boundary, the 
increase may be suppressed. However, this is a very difficult back-reaction problem. The boundary conditions are normally 
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used to determine the inner structures. The adjustable boundary conditions should be controlled by the interior. Thus, the 
boundary conditions and the inner structures should be determined simultaneously. Such numerical scheme is not known and 
should be developed in future. Otherwise, extensive study to find out consistent boundary conditions for large rj is required. 
Our numerical method will be improved by either approach. 

Pulsar magnetosphere is described in most spatial region by the force-free and ideal MHD approximations, which corre- 
spond to the limit of 77 ^ 1. Stationa ry axisymmetric magnetosphere is constructed so far by a solutio n of the Grad-Shafranov 
equation with these approximations ( Contopoulos et al. 19991 : Oeura fc Koiimal 20031 : Gruzinov 2005 ) . Poloidal magnetic field 
approaches a quasi-spherical wind at infinity. Ther e is a discontin u ity in toroidal magnetic field at the boundary of the coro- 
tating region, where the current sheet is formed. iLovelace et al. l|2006t ) have obtained an alternative solution of the same 
equation, but with a different injection current. Their solution exhibits a jet along polar axis and a disk on the equator. 
There is no current sheet in their numerical model. Thus there are at least two models in the strong magnetic field limit, 
?7 S> 1. The poloidal magnetic field at infinity is quasi-spherical and there is no current sheet in our numerical solution. It is 
interesting to examine the model sequence by increasing the parameter rj. Some plasma flow should be highly constrained to 
a thin region, or the jet-disk system should be formed in the large rj limit. It is unclear whether or not many solutions exist, 
depending on physical situations including the plasma state. In order to address these questions, it is necessary to construct 
the magnetosphere with plasma flow beyond the force-free and ideal MHD approximations. We have here presented a possible 
approach, although the improvement is needed. 
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APPENDIX A: DIMENSIONLESS FORMS 

We consider the dimensionless forms of eqs. (|3]),([5l and (|18p . The magnetic function G is normalized in terms of the magnetic 
dipole moment ^ and the distance to light cylinder tl as G = fiG^ /rL, where a symbol ^ denotes a dimensionless quantity. 
The electric potential is expressed as = /if2^'^/(cr_L) = ji^'^ /r\. Two integrals (|15|) and (|16p are written as J± — 7ocr_LjJ^ 
and K± — ^oK\, where 



7± 
70 

7± 
70 
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(Ai) 

(A2) 



These values depend on a dimensionless parameter t] = q^/{m(?^or\). As for the stream function F±, the normahzation con- 
stant is /if2/(qrz,) = fic/{qr\) and F± = ^,Q,F]^/ {qrL)- The number density is normalized from ea. (|12|) as n± — ^VLn^_^ / {qcr\) . 



Electric charge and current densities are normalized as pe = 
functions, eqs. Q,© and (|18|) can be written as 

V^ln 



^Q,p],/(cr\) and j = /ifij^/rl, . Using these dimensionless 

(A3) 
(A4) 



n^7o 
7± 



+ 



c^7o 



7±K.+«±e)V ± ci?t i 



7± 



(A5) 



where i?^ — R/tl, r'' — r/rL and the differential operators with the symbol ^ are defined by — r/rL- From these expressions, 
we find that rj is an important parameter. 



